%% Markov Chain Monte Carlo (MCMC)

N=10000;  %number of steps in  chain

%Intitialize chain vector
Chain=zeros(N,1);

%Setup initial value
Chain(1)=10;

%Setup probability function to explore
Mu=0;
Sigma=1;
p_x=@(x)normpdf(x,Mu,Sigma);

%Setup test jump size
TestSigma=.1;

%run chain

for nn=1:N-1
    %Make test jump
    CurrentX=Chain(nn);
    TestX=CurrentX+TestSigma*randn();
    
    %Metropolis-Hastings rule
    if rand()< (p_x(TestX)/p_x(CurrentX))
        Chain(nn+1)=TestX;
    else
        Chain(nn+1)=CurrentX;
    end
end

figure
hist(Chain,20)
xlabel('x')
ylabel('p(x)')

%% Visualize Chain

Xleft=-3;
Xright=11;

Xvec=(Xleft:.1:Xright);
Yvec=1:N;
FigH=figure

for nn=N/2:N
    figure(FigH)
    %plot chain
    subplot(2,1,1)
    plot(Chain(1:nn),Yvec(1:nn))
    hold on
    plot(Chain(nn),Yvec(nn),'ro','Linewidth',3)
    axis([Xleft Xright 0 N])
    hold off
    xlabel('x')
    ylabel('Chain Index')
    
    %plot histogram
    subplot(2,1,2)
    [Counts]=hist(Chain(1:nn),Xvec);
    bar(Xvec,Counts)
    axis([Xleft Xright 0 max(Counts)*1.1])
    xlabel('x')
    ylabel('p(x)')
    
    pause(.001)
    
    
end



